Recoil velocities from equal-mass binary black-hole mergers: a systematic investigation of 

spin-orbit aligned configurations 
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Binary black-hole systems with spins aligned with the orbital angular momentum are of special interest, as 
studies indicate that this configuration is preferred in nature due to non-vacuum environmental interactions, 
as well as post-Newtonian (PN) spin-orbit couplings. If the spins of the two bodies differ, there can be a 
prominent beaming of the gravitational radiation during the late plunge, causing a recoil of the final merged 
black hole. In this paper we perform an accurate and systematic study of recoil velocities from a sequence 
of equal-mass black holes whose spins are aligned with the orbital angular momentum, and whose individual 
spins range from a — +0.584 to —0.584. In this way we extend and refine the results of a previous study 
which concentrated on the anti-aligned portion of this sequence, to arrive at a consistent maximum recoil of 
448 ± 5 km/s for anti-aligned models as well as to a phenomenological expression for the recoil velocity as 
a function of spin ratio. Quite surprisingly, this relation highlights a nonlinear behavior, not predicted by the 
PN estimates, and can be readily employed in astrophysical studies on the evolution of binary black holes in 
massive galaxies. An essential result of our analysis, without which no systematic behavior can be found, is the 
identification of different stages in the waveform, including a transient due to lack of an initial linear momentum 
in the initial data. Furthermore, by decomposing the recoil computation into coupled modes, we are able to 
identify a pair of terms which are largely responsible for the kick, indicating that an accurate computation can 
be obtained from modes up to ^ = 3. Finally, we provide accurate measures of the radiated energy and angular 
momentum, finding these to increase linearly with the spin ratio, and derive simple expressions for the final spin 
and the radiated angular momentum which can be easily implemented in A^-body simulations of compact stellar 
systems. Our code is calibrated with strict convergence tests and we verify the correctness of our measurements 
by using multiple independent methods whenever possible. 

PACS numbers: 04.25.Dm, 04.30.Db, 95.30.Sf, 97.60.Lf 



I. INTRODUCTION 

Recent developments in numerical relativity have solved 
the problem of stably evolving black hole initial data for use- 
ful timescales, and opened the door to studies of physical phe- 
nomena resulting from strong-field gravitational interactions. 
A result of particular interest to astrophysics is an accurate 
calculation of the recoil velocity which is generated during an 
asymmetric collision of a black-hole binary. It is well known 
that a binary with unequal masses or spins of the individual 
bodies will radiate gravitational energy asymmetrically. This 
results in an uneven flux, which gives a net linear momentum 
to the final black hole, often called a "kick" 01 0]. While 
estimations of kick velocities have been available for some 
time iHHIlt], the largest part of the system's acceleration is 
generated in the final orbits of the binary system, and as such 
requires fully relativistic calculations to be determined accu- 
rately. 

Over the past year, a number of numerical relativity simu- 
lations have been carried out to determine recoil velocities in 
various sections of the parameter space of binary black-hole 
systems. The first systems to be studied were unequal mass 
systems with moderate mass ratios, where the first calcula- 
tions were performed by the Penn State [@] and Goddard 1 7| 
groups, with simulations at mass ratios near the estimated 



peak of the Fitchett formula |3l ■ A more extensive study, ex- 
ploring a large number of models between mass ratios 0.25 
to 1.0, was carried out by the Jena group iH, providing for 
the first time a mapping of the unequal mass parameter space 
with fully relativistic simulations. The recoils from systems in 
which the bodies had spin were first considered by a number 
of studies in the first half of this year. The Penn State group 
examined a sequence of equal mass binaries with spins equal 
and anti-aligned, determining that the largest recoilpossible 
from such an evolution is of the order of 475 km/s At the 
same time, in Ref. ifioll we studied a sequence of models in 
which the spins are anti-aligned, but of different magnitude, 
and arrived at a similar estimate of 450 km/s. The Jena and 
Brownsville (now Rochester) groups showed that extremely 
large kicks are possible from particular configurations of mis- 
aligned spins, measuring recoils as high as 25001<;m/s ifTTIl . 
and extrapolating to 4000 km/s for the maximally spinning 
case idlil. Such spin configurations have recently been 
studied in more detail in fl4j]. Velocities of this magnitude 
have a number of astrophysical implications for models of 
galaxy mergers. 

In this paper, we expand on the work performed in ifioll . 
extending it in a number of different ways. First, we con- 
sider a larger sequence of aligned but unequal spins with spin- 
ratios ranging from —1 to +1, where the spins are aligned (or 
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anti-aligned) with the orbital angular momentum. Our inter- 
est in this set of models is motivated by the fact that there are 
strong indications that binary black-hole systems having spins 
aligned with the orbital angular momentum are preferred in 
nature. Post-Newtonian (PN) studies in vacuum have in fact 
shown that in vacuum the gravitational spin-orbit coupling has 
a tendency to bring about such an alignment from generic ini- 
tial conditions [15]. Furthermore, in astrophysical situations 
where there is likely to be at least some component of inter- 
stellar matter inducing a dissipative dynamics, there is also a 
tendency to align lfl6ll . 

We describe the influence of the initial dynamics on the ra- 
diated waveforms and the importance of suitable vector inte- 
gration constants to remove these effects when determining 
the final recoil velocity. These vectors, in fact, capture the in- 
formation about the net linear momentum that the spacetime 
has built-up during its past evolution and prior to the actual nu- 
merical evolution and can result into a significant correction. 
We discuss how to use the results obtained to derive a phe- 
nomenological expression for the recoil velocity as a function 
of the spin ratio. Finally, we also compute how the angular 
momentum of the system is redistributed between radiation 
and the spin of the final black hole, providing useful expres- 
sions as functions of the spin ratio. 

The paper is organized as follows: Section|II]describes the 
code, as well as the initial data construction, and calibration 
tests. In Sectionlllllwe discuss the calculation of the recoil ve- 
locity from gravitational-wave data on a large sphere. We in- 
troduce and compare two methods, one based on the Newman- 
Penrose ^'4 scalar which is the usual method that has been 
adopted in recent numerical studies, and another which is 
based on perturbations of Schwarzschild black holes modeled 
by a gauge-invariant formalism. Though the two methods are 
based on quite different underlying assumptions, they agree 
very well in their estimation of physical quantities, and in par- 
ticular the recoil velocity. Section |IV] describes evolutions 
of the ahgned-spin sequence and the dependence of the re- 
coil velocity on the spin-ratio. We find that the data show 
an almost linear behavior at large negative spin-ratios, as pre- 
dicted by PN calculations. However taking into account also 
results from positive spin-ratios, the data suggest a nonlin- 
ear (quadratic) dependence and we give a phenomenological 
expression for the recoil velocity as a function of the spin ra- 
tio. Extrapolating our results to the case of maximally rotat- 
ing black holes, we find that the maximum recoil velocity at- 
tainable by spin-orbit aligned configurations is 448 ± 5 km/s. 
Finally, we discuss the radiation of mass and angular momen- 
tum for these evolutions, determining the parameters of the 
isolated final black holes and show the excellent conservation 
of mass and angular momentum recorded in our simulations. 
Again we provide phenomenological expressions for the rel- 
ative amount of radiated mass and spin as functions of the 
initial spin ratio. 

In the following equations we use Greek indices (running 
from to 3) to denote components of four-dimensional objects 
and Latin indices (running from 1 to 3) for three-dimensional 
ones that are defined on space-like foliations of the space- 
time. 



II. MATHEMATICAL AND NUMERICAL SETUP 

The data presented in this paper were produced using the 
CCATIE code, a three-dimensional finite differencing code 
based on the Cactus Computational Toolkit iflTiflsll . The cur- 
rent code is an evolution of previous versions which imple- 
mented an excision method and co-rotating coordinates lfl9l 
H^IIH]. The main features of the code, in particular the evo- 
lution equations, remain the same. However, some modifi- 
cations have been introduced in the gauge evolution to ac- 
commodate "moving punctures" which has proven to be an 
effective way to evolve black hole spacetimes |22, 23]. This 
method simply removes any restrictions on movement of the 
punctures from their initial locations, allowing them to be ad- 
vected on the grid. 



A. Evolution system 

We evolve a conformal-traceless "3 + 1" formulation of 
the Einstein equations lfl9l l24l l25i l26ll . in which the space- 
time is decomposed into three-dimensional spacelike slices, 
described by a metric 7^ , its embedding in the full spacetime, 
specified by the extrinsic curvature Kij, and the gauge func- 
tions a (lapse) and /3' (shift) that specify a coordinate frame 
(see Sect. Ill Bj for details on how we treat gauges and IztIi for 
a general description of the 3 + 1 split). The particular system 
which we evolve transforms the standard ADM variables as 
follows. The 3-metric 7^ is conformally transformed via 



= TIT lndet7i-,-, 7y = e (1) 



and the conformal factor (f> evolved as an independent vari- 
able, whereas 7^ is subject to the constraint det 7^ = 1. The 
extrinsic curvature is subjected to the same conformal trans- 
formation, and its trace tr Kij evolved as an independent vari- 
able. That is, in place of Kij we evolve: 

K = tr lu, = g'^ lUj , = (Xy - ^7,, K) , (2) 

with tr Aij — 0. Finally, new evolution variables 

r = 7^"f (3) 

are introduced, defined in terms of the Christoffel symbols of 
the conformal 3-metric. 

The Einstein equations specify a well known set of evolu- 
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tion equations for the listed variables and are given by 



{dt- Cp)%j=-2aA,j, (4) 
{dt~ C0)<i>^~^aK, (5) 

+ -2i,fci^), (6) 



{dt - Cp) K = -DW^a + a{A,jA'^ + -A'^), (7) 

- 2A'^dja + 2a(Pjfci^'= + GA'^djcj) 

- IrdjK), (8) 

where Rij is the three-dimensional Ricci tensor, Di the co- 
variant derivative associated with the three metric jij and 
"TF" indicates the trace-free part of tensor objects. The Ein- 
stein equations also lead to a set of physical constraint equa- 
tions that are satisfied within each spacelike slice, 

n EE ^ _ K,jK'^ ^ 0, (9) 

M' = Dj {K'^ - f^K) = 0, (10) 

which are usually referred to as Hamiltonian and momentum 
constraints. Here i?'^' = RijY^ is the Ricci scalar on a three- 
dimensional time slice. Our specific choice of evolution vari- 
ables introduces five additional constraints, 

dct7,j = l, (11) 
tri,.,=0, (12) 

r-7^*^fj,. (13) 

Our code actively enforces the algebraic constraints (fTTl i 
and (fT2] i. The remaining constraints, H, M.\ and ( flJl l, are not 
actively enforced, and can be used as monitors of the accuracy 
of our numerical solution. See || 20il for a more comprehensive 
discussion of the these points. 

B. Gauges 

We specify the gauge in terms of the standard ADM lapse 
function, a, and shift vector, /3° ll28ll . We evolve the lapse 
according to the "1 + log" slicing condition: 

dta- P'd,a = -2a{K - Ko), (14) 

where Kq is the initial value of the trace of the extrinsic cur- 
vature, and equals zero for the maximally sliced initial data 
we consider here. The shift is evolved using the hyperbolic 
f -driver condition ll20ll . 

dtP'-P'djl3' = ^aB\ (15) 
dtB'-P^djB' = dtT' - P^ djT' - ijB\ (16) 



where 77 is a parameter which acts as a damping coefficient. 
The advection terms on the right-hand-sides of these equa- 
tions were not present in the original definitions of lEoll . where 
co-moving coordinates were used, but have been added fol- 
lowing the experience of ll29i[30ll . and are required for correct 
advection of the puncture in "moving-puncture" evolutions. 



C. Numerical methods 

Spatial differentiation of the evolution variables is per- 
formed via straightforward finite-differencing using fourth- 
order accurate centered stencils for all but the advection terms 
for each variable, which are upwinded in the direction of the 
shift. Vertex-centered adaptive mesh-refinement (AMR) is 
employed using nested grids ll3ll [32I1 with a 2 : 1 refinement 
for successive grid levels, and the highest resolution concen- 
trated in the neighborhood of the individual horizons. Individ- 
ual apparent horizons are located every few time steps during 
the evolution 

The time steps on each grid are set by the Courant condi- 
tion and thus the spatial grid resolution for that level, with the 
time evolution being carried out using fourth-order accurate 
Runge-Kutta integration steps. Boundary data for finer grids 
are calculated with spatial prolongation operators employing 
5th-order polynomials, and prolongation in time employing 
2nd-order polynomials. The latter allows a significant mem- 
ory saving, requiring only three time levels to be stored, with 
little loss of accuracy due to the long dynamical timescale rel- 
ative to the typical grid time step. 

In the results presented below we have used 8 levels of 
mesh refinement with finest grid resolutions of h/M ~ 0.030, 
0.024, and 0.018; we will refer to these resolutions as "low", 
"medium" and "high" respectively. We find that the medium 
(i.e., h = 0.024 M) fine-grid resolution is typically good 
enough to accurately represent the dynamics which we are 
studying here and will be used hereafter as our fiducial res- 
olution. In this case, the wave-zone grid has a resolution of 
h = 1.536 A/. In addition, when measuring the convergence 
order (see discussion in Sect. Ill Ek we have also used a "very- 
high" resolution of h/M — 0.012 which therefore gives a 
factor of 2 refinement with respect to the "medium" resolu- 
tion; this should be contrasted with similar convergence tests 
recently discussed in the literature and in which the refinement 
factor is much smaller. 

The finest grids are centered on each black hole, with a 
radius about 50% larger than the apparent horizon. A sin- 
gle grid resolution covers the region between r = 20 M and 
r = 80 M, in which our wave extraction is carried out. The 
outer (coarsest) grid extends to a spatial position which is 
large compared with the evolution time of the system. In 
particular, it ranges from 256 M in each coordinate direction 
for the binaries which merge rapidly, up to 768 M for the bi- 
naries which inspiral more slowly because of the spin-orbit 
interaction. In all cases, artificial wave-like boundary condi- 
tions are used, and although these are not explicitly constraint- 
preserving, they do not introduce major violations of the con- 
straints as long as they are placed sufficiently far away from 
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the central black holes {i.e., with a light-crossing time which 
is large as compared to the time for the merger). Further- 
more, for the models considered here, in which all spins are 
directed along the z-axis of our Cartesian grids, it is possible 
to use a reflection symmetry condition across the z = plane. 
Tests against the runs on a full grid show that this symmetry 
is preserved to a high degree in our simulations (i.e., with dif- 
ferences below 10^^"*) so that this symmetry boundary has no 
influence on the dynamics. 

D. Initial data 

The initial data are constructed applying the "puncture" 
method Issll . which uses Bowen-York extrinsic curvature and 
solves the Hamiltonian constraint equation numerically as 
in ill. 

We have considered a sequence of binaries for which the 
initial spin of one of the black holes is held fixed at ^2 /^/'^ = 
0.146 e^, and the spin of the other black hole is Si/Al'^ = 
(ai/a2)S2/-A/^, where the spin ratio ai/a2 takes the values 
-1, -3/4, . . . , 3/4, 1, and M is the sum of the black hole 
masses, M = Mi + Al2. Thus the black hole spins are anti- 
aligned when 01/02 is negative and aligned when it is posi- 
tive. In all cases the initial data parameters are chosen such 
that the black hole masses are 

Fa, AnSf 1 

IstIIssIi where Ai is the area of the i-th apparent horizon. 

For the orbital initial data parameters we use the effective 
potential method introduced in ll39ll and extended to spinning 
configurations in 1I40I1 . The effective potential method is a way 
of choosing the initial data parameters such that the required 
physical parameters (e.g. masses and spins) are obtained to 
describe a binary black-hole system on a quasi-circular orbit. 

The free parameters to be chosen for the puncture initial 
data are: the puncture coordinate locations Ci, the puncture 
mass parameters nii, the linear momenta p^, and the individ- 
ual spins Sj . Since we are interested in quasi-circular orbits 
we work in the zero momentum frame and choose = — 
to be orthogonal to C2 — Ci. The physical parameters we 
want to control are: the blackhole mass ratio M1/M2, the or- 
bital angular momentum L = Ci x p^ + C2 ^ P2 (see for 
example iHliSlilll) and the dimensionless spin parameters 
Ci = Si/Mf. In order to choose the input parameters that 
correspond to the desired physical parameters we have to use 
a non-linear root finding procedure, since the physical param- 
eters depend non-linearly on the input parameters and it is not 
possible to invert the problem analytically. 

As detailed in ll40ll . when the black-hole spins are taken 
as parameters, it is possible to reduce the number of inde- 
pendent input variables, so that at a given separation C = 
IC2 — Ci|/toi, the independent input parameters are: q = 
mi/ 1712 and the dimensionless magnitude of the linear mo- 
mentum p/mi. Using a Newton-Raphson method, we solve 
for q and p/mi so that M1/AI2 = 1 and the system has 



a given dimensionless orbital angular momentum, L/{^M) 
where fj, = mim2/Af^ is the reduced mass. For such a con- 
figuration the initial data solver ll36ll returns a very accurate 
value for Af^^^j, which together with the accurate irreducible 
mass calculated by the apparent horizon finder ll33ll42ll makes 
it possible to calculate an accurate value of the dimensionless 
binding energy 

Et,/^l = (A/^„„ - Ml - M2)/y.. (18) 

The quasi-circular initial data parameters are then obtained by 
finding the minimum in Ei,/ fjL for varying values of C while 
keeping the required orbital angular momentum L/{iiM) 
constant. 

We chose a fixed orbital angular momentum L/{fiAI) = 
3.3 for our quasi-circular orbit initial data parameters. This 
value was chosen to ensure that model rO would have enough 
evolution time for an accurate kick measurement, while at the 
same time model r8 would not require too much evolution 
time. In order to check the influence of the evolution time 
before plunge on the kick measurements of the rO model, we 
also calculated initial data for a rO configuration at larger ini- 
tial separation rOl and at smaller initial separation rOs. The 
parameters for all the initial data sets are shown in Table |T] 

Note that the physical mass Mi of a single puncture black 
hole increases when the spin parameter is increased if the 
mass parameter rrii is kept constant. For that reason obtain- 
ing Ml = M2 in general requires that mi 7^ 7TI2. Even in 
the case where the spins have the same magnitude but differ- 
ent directions, the two black holes will have different spin- 
orbit interactions leading to slightly different physical masses 
if mi = m2. For this reason, the initial data for rO in Table U 
has slightly different puncture mass parameters mi 7^ m2- 
In contrast, in model r8 the black holes have identical spin 
parameters and thus also the same spin-orbit interaction, re- 
sulting in identical mass parameters mi = m2. 

E. Convergence tests 

As described in Section III CI the finite difference error 
of the derivative stencils used in the numerical algorithm is 
0{h'*'), while the error in the time-interpolation stencils used 
for mesh refinement boundary points is 0{At^). Thus the 
expected theoretical convergence rate is three. However, it 
is only time-related operations which are at third order, and 
since the time step which we use is smaller than the grid 
spacing and much smaller than the dynamical timescales, we 
can expect that the error coefficient of the leading order term 
is quite small. Third order convergence is expected during 
time-periods when the system goes through rapid dynamical 
changes, such as the plunge or merger. 

The proper convergence of the code was established using 
the binary system rO, for which we have carried out evolutions 
using 8 levels of mesh refinement with fine grid-spacings of 
h/M = 0.024, 0.018, and 0.012 (i.e., resolutions "medium", 
"high", and "very-high", respectively, where "low" refers to 
h — 0.030 which was deemed to be of insufficient accuracy 
for the results of this paper). Other refinement levels have 
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TABLE I: The puncture initial data parameters defining the binaries: location ±x/M, linear momenta ±p/M, mass parameters rrii /A/, spins 
Si/M"^ , dimensionless spins ai, ADM mass M^^, measured at infinity, and ADM angular momentum J^dm computed from Eq. J47t . Note 
that we set Mi = M2 = 1/2 [cf., Eq. iITtIiI. 



Model 


±x/A/ 


±p/M 


mi /M 


7712 /M 


Si/M^ 


S2/M^ 


ai 


12 






rO 


3.0205 


0.1366 


0.4011 


0.4009 


-0.1460 


0.1460 


-0.5840 


0.5840 


0.9856 


0.8252 


rl 


3.1264 


0.1319 


0.4380 


0.4016 


-0.1095 


0.1460 


-0.4380 


0.5840 


0.9855 


0.8612 


r2 


3.2198 


0.1281 


0.4615 


0.4022 


-0.0730 


0.1460 


-0.2920 


0.5840 


0.9856 


0.8979 


r3 


3.3190 


0.1243 


0.4749 


0.4028 


-0.0365 


0.1460 


-0.1460 


0.5840 


0.9857 


0.9346 


r4 


3.4100 


0.1210 


0.4796 


0.4034 


0.0000 


0.1460 


0.0000 


0.5840 


0.9859 


0.9712 


r5 


3.5063 


0.1176 


0.4761 


0.4040 


0.0365 


0.1460 


0.1460 


0.5840 


0.9862 


1.007 


r6 


3.5988 


0.1146 


0.4638 


0.4044 


0.0730 


0.1460 


0.2920 


0.5840 


0.9864 


1.044 


r7 


3.6841 


0.1120 


0.4412 


0.4048 


0.1095 


0.1460 


0.4380 


0.5840 


0.9867 


1.081 


r8 


3.7705 


0.1094 


0.4052 


0.4052 


0.1460 


0.1460 


0.5840 


0.5840 


0.9872 


1.117 


rOZ 

rOs 


4.1924 
2.8186 


0.1073 
0.1441 


0.4066 
0.3997 


0.4065 
0.3994 


-0.1460 
-0.1460 


0.1460 
0.1460 


-0.5840 
-0.5840 


0.5840 
0.5840 


0.9889 
0.9849 


0.8997 
0.8123 



resolutions that are half of the next finest grid. The refinement 
levels on the initial slice are set up to be identical for the three 
resolutions and their locations and sizes evolve according to 
the same algorithm in each case. 

We focus on the convergence of a number of different as- 
pects of the code. The first of these is the degree of satis- 
faction of the Einstein equations, which can be partially de- 
termined by examining the Hamiltonian and momentum con- 
straints (l9])-(fTQli. A more stringent requirement is to evaluate 
how well the Einstein tensor satisfies the vacuum condition, 
GajS = 0. For this we define the positive definite quantity 

Q — \ V + ^01 + • • • + ^"§3 outside appar horizons 
~ |_ inside appar horizons . 

(19) 

In computing norms over the entire grid, we find it useful to 
mask out the interiors of the horizons, where the error at the 
puncture locations - which is not expected to converge - can 
dominate over more relevant errors in the physically observ- 
able domain. In order to compute Gap we compute the 4- 
derivatives of the ADM metric, lapse and shift, then construct 
the 4-derivatives of the 4-metric from which we can com- 
pute the Riemann tensor and then finally obtain Gap- Time- 
derivatives are taken using three time-levels, centered around 
the past time-level. Spatial derivatives are taken using fourth- 
order accurate centered stencils. Thus the finite-difference er- 
ror in computing Gap is 0{AP) in time and 0{h'^) in the 
space dimensions. Effectively we see a minimum of third or- 
der accuracy for this quantity, indicating that the coefficient of 
the 0{At^) error term is small compared to the higher-order 
terms. 

Since the metric gradients and hence the truncation errors 
are the largest near the black-holes, through the Loo norm 
of iT% we effectively monitor that the Einstein tensor con- 
verges near the horizons for the duration of the evolution. We 
regard this as a rather stringent test in comparison with the 
common use of the L2 norm, as the latter tends to dilute er- 
rors in small regions or 2D surfaces such as grid boundaries, 
as they are normalized over the entire grid volume. By con- 
trast, the Loo norm measures the worst error on the grid, which 
by propagation of error will also suffer if there are any non- 
convergent regions on the grid. 



This convergence of G is summarized in Fig. [T] which re- 
ports the time evolution of the Loo norm of iT% at the medium 
and very-high resolutions. Also indicated with dashed and 
dotted lines are the expression for the Loo norm of ( fT9] l at the 
very-high resolution when rescaled for third (dotted line) and 
fourth-order convergence (dashed line). 

There is a period at the beginning of the evolutions where 
the initial data construction prevents fourth-order conver- 
gence. This is due to the fact that the initial data is computed 
by an interpolation of the results of a spectral solver onto the 
finite difference grid which is used for evolution. An error 
is introduced because we keep fixed the number of spectral 
coefficients and because the Cartesian grid points do not co- 
incide with the spectral collocation points of the Chebyshev 
polynomials, resulting in a certain amount of high-frequency 
noise that spoils the convergence for some time at the begin- 
ning of the simulation. Numerical dissipation and the con- 
straint damping built into the evolution system implies that 
the evolution quickly adjusts itself to actually solving the Ein- 
stein equations to a good accuracy. The effects of these initial 
transient modes can last for different amounts of time for the 
different resolutions, e.g., ^ 10 M for the medium resolution 
and ~ 30 AI for the very-high resolution. 

Soon after this transient has disappeared, the code shows 
the expected fourth-order convergence, with the largest values 
of the violation found in the vicinity of the apparent horizons, 
where the gradients in the metric are the steepest. The vio- 
lations grow rapidly with time as the binary inspirals and the 
largest values of the violation of the Einstein tensor are seen 
at the time of the merger, i « 109 AI, with values as large as 
0(300). Such violations are essentially confined to a single 
grid point on the trailing edge of the apparent horizon and are 
produced by the very steep gradients in the shift. Clearly, vi- 
olations of this magnitude would not be revealed when look- 
ing at the L2 norms and are a source of concern. However, 
as we will show later, such violations do not propagate away 
from the horizon to affect the fourth-order convergence of the 
waveforms. 

At the time of the merger the excision of a common appar- 
ent horizon from the calculation of the Loc norm is respon- 
sible for the decrease by about four orders of the violation. 
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FIG. 1: The Loo norm of the Einstein tensor Eq. l |19t as a func- 
tion of time. During the periods of strong dynamics (i.e., when 
the time derivatives of the evolution variables are large) the conver- 
gence order is dominated by the accuracy of the time-interpolation 
algorithm used at mesh refinement boundaries, thus yielding third- 
order accuracy. At the times when these time-derivatives are small, 
the fourth-order finite-differencing algorithm becomes the dominant 
source of the error. Note that the very large violations (of 0(300) 
at the medium resolution) are confined to a single grid point on the 
trailing edge of the apparent horizon and are produced by the very 
steep gradients in the shift. As discussed later, this does not affect 
the fourth-order convergence of the waveforms. At the time of the 
merger a common apparent horizon forms and its excision from the 
calculation of the Loo norm is responsible for the drop in the viola- 
tion. 



After this, the Loo do not grow further in time for the very- 
high resolution simulation, while a modest increase is seen in 
the simulation run at medium resolution. During this time the 
code shows a convergence which is between third-order (right 
after the merger) and fourth-order (during the ringdown). 

In addition to convergence in the Einstein tensor, we also 
validate the correctness of the physically relevant information 
contained in the waveforms. We do this by computing con- 
vergence rate of the waveforms Q22, Qtd^ ™d Q21 using the 
ratio of the integrated differences between the medium and 
high resolutions, and the high and very-high resolutions 



024 ■ 



^0.0181 



0.018 



Qo.012prf" 



(20) 



where u = i — r^, is the retarded time at a given detector, 
Q stands for either Q^j' Q33 or Q21 ™d refers to either its 



amplitude or the phase. As indicated in Eq. (l20b . the integrals 
are evaluated over the retarded interval [wi , U2] which does not 
include the initial spurious burst of radiation (which we do 
not expect to converge) but contains otherwise the complete 
waveform including the ringdown. 

Assuming a truncation error 0{hP) and that the coefficient 
of this error does not depend on resolution, the function p be- 
comes to leading order 



(^10.024)^ — (/^0.018)^ 

{ho.oisy - {ho.oi2)P 



(21) 



where /io.o24 — 0-024 M and we we underline the importance 
of having used a full doubling of the resolution between the 
smallest and largest resolution to improve the accuracy of this 
estimate over more narrowly spaced resolution steps. In prac- 
tice, we measure p and then solve for the "effective " conver- 
gence order p using Eq. ( I2TI ). A discussion of the details in 
this procedure are presented in Appendix|A]alongside with the 
computed convergence rates for the amplitudes and phases of 
Q which are found to between two {£ = 3) and four {£ = 2) 
(cf., Tableig. 

It should be noted that the above definition of convergence 
rate naturally results in non-integer values for the exponent p, 
even though our methods are explicitly polynomial. This is 
because the derivation of ( 1211 1 assumes a coefficient of one in 
the leading order error term that extrapolates between the res- 
olutions. If the coefficient is in practice different for a given 
set of resolutions, then a non-integer value results which is 
larger if the coefficient is smaller. As such, values obtained in 
this way should not be considered literal polynomial extrap- 
olation orders. By "convergence order 3.8" we rather mean 
that our results are consistent with third-order finite differ- 
encing where the leading third-order error coefficient is quite 
small so that at the given resolutions the convergence appears 
to be closer to a fourth-order approximation. Very high con- 
vergence exponents are a likely indication that the lowest res- 
olution is not in the convergent regime for the measured quan- 
tity. Non-integer convergence orders obtained in this way are 
resolution dependent, and should themselves converge to the 
lowest order finite difference approximation used in the code 
in the limit of infinite resolution. 

An important property of the waveforms which has 
emerged when performing these convergence tests is that the 
dominant source of error is a de-phasing which causes the 
lower resolution evolutions to "lag" behind the higher reso- 
lution. This delay is usually rather small and between 0.1 M 
and 0.5 M, but it is clearly visible when comparing the total 
amplitude of Q as a function of time. The most important 
consequence of this error is that it can spoil the convergence 
tests if not properly taken into account: the residuals errors 
seem, in fact, to indicate over-convergence. This is shown in 
the upper panel of Fig. |2l which reports the differences be- 
tween Q22 when computed at different resolutions scaled for 
fourth-order convergence. Clearly the overlap is rather poor 
and even indicating that the truncation error is smaller than 
expected. This is obviously an artifact of the near cancellation 
of the lowest-order terms in the truncation error and induced 
by the small time-differences at different resolutions. 
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FIG. 2: Convergence of the fiducial waveform Q for binary sys- 
tem rO before and after tfie time-shift defined in Eqs. ( IAlb -( lA3l l. In 
the upper graph we show the difference between Q22 when computed 
at different resolutions, scaled for fourth-order convergence and us- 
ing raw data (i.e., without time-shifting). The overlap between the 
curves is rather poor indicating an over-convergence (i.e., the trunca- 
tion error appears to be smaller than expected). In the lower panel we 
show the same data but after time-shifting. The very good overlap of 
the scaled curves on the indicates that the time-shifting is essential 
for obtaining properly scaling differences between runs of various 
resolutions. 



FIG. 3: Accuracy of the fiducial waveform Q22 for the binary sys- 
tem rO. In the upper graph we show the waveforms at the three dif- 
ferent resolutions: very-high (continuous line), high (dashed line), 
medium (dotted line). The accuracy is very good already with the 
lowest resolution and the curves cannot be distinguished. The lower 
panels show magnifications of some relevant portions of the wave- 
form, with the lower-left panel concentrating on the initial transient 
radiation produced by the truncation error. The lower-right panel, on 
the other hand, refers to the quasi-normal ringing and shows that it is 
well-captured at all resolutions. 



We remove this effect by shifting the time coordinate of the 
medium and high resolution runs by the time interval needed 
to produce an alignment of the maxima of the emitted radia- 
tion. Details on how to do this are discussed in Appendix lAl 
and we report in the lower panel of Fig. |2] the same data 
shown in the upper panel, but after the time-shifting. Clearly, 
the overlap is now extremely good suggesting that the time- 
shifting is essential for obtaining the expected fourth-order 
convergence in the waveforms. In accord with the conver- 
gence in the waveforms we also see fourth order convergence 
in the final kick value. 

As a final note we remark that besides validating a proper 
convergence of the code, it is also important to assess the ac- 
curacy of any measurable quantity at the relevant resolutions 
considered here. As a representative and physically mean- 
ingful quantity we have considered the accuracy of the fidu- 
cial waveform Q22 for the binary system rO. This is shown 
in Fig. [3] where in the upper graph we report the waveforms 
at the three different resolutions: very-high (continuous line), 
high (dashed line) and medium (dotted line). Already with 
the lowest of these resolutions the accuracy is sufficiently high 



so that the curves are essentially indistinguishable from each 
other by eye. The lower panels show magnifications of the 
relevant portions of the waveform, with the lower-left panel 
concentrating on the initial transient radiation produced by 
the truncation error. The latter clearly is rather large at the 
medium resolution, but it nicely converges away when the 
grid spacing is decreased. The lower-right panel, on the other 
hand, refers to the quasi-normal ringing and shows that it is 
well-captured at all resolutions. 



III. LINEAR MOMENTUM OF BLACK HOLE 
SPACETIMES 

In radiating spacetimes where the radiation is emitted 
asymmetrically, there will be a net linear momentum imparted 
to the system. In particular, in the case of a binary black hole 
merger, the final black hole receives a "kick" which causes 
it to move off at a given velocity. This velocity can be de- 
termined by an analysis of the emitted radiation. In ADM- 
type numerical simulations, this is typically done by evalu- 
ating some scalar quantity which can be associated with the 
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wave energy at some large radius within the computational 
domain. The chosen radius needs to be large enough that it 
is in the "wave zone", where non-linear self-interaction of the 
gravitational field is negligible and the waves can be picked 
out as perturbations of a background. 

Two methods have become commonplace to determine the 
emitted wave energy. The first uses the Newman-Penrose cur- 
vature scalar 4*4, which can be identified with the gravitational 
radiation if a suitable frame is chosen at the extraction radius. 
An alternative method measures the metric of the numerically 
generated spacetime against a fixed background at the extrac- 
tion radius, and determines the Zerilli-Moncrief perturbation 
modes. Both methods yield data for the gravitational wave 
energy which can be integrated to determine a net linear mo- 
mentum, as described in more detail in the following sections. 



A. Kick measurements via ^'4 

The Newman-Penrose formalism provides a convenient 
representation for a number of radiation related quantities as 
spin-weighted scalars. In particular, the curvature component 
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^4 



(22) 



is defined as a particular component of the Weyl curvature, 
CapfS, projected onto a given null frame, {I, n, m, fh}. In 
practice, we define an orthonormal basis in the three space 
{r,6,4)), centered on the Cartesian grid center and oriented 
with poles along z. The normal to the slice defines a time-like 
vector t, from which we construct the null frame 



FIG. 4: Amplitude of r^, ^^j^ |^'4| for extraction spheres at = 
30 M, 40 M, 50 M and 60 M, demonstrating that *4 does indeed 
fall off as required by the peeling property. There is a slight de- 
crease in amplitude with larger radius, suggesting that dissipative ef- 
fects may become important at larger radii. Results in this paper use 
wavefoiTns from the — 50 AI extraction sphere, unless indicated 
otherwise. 



m= -^{6 -let)) 



(23) 

We then calculate ^'4 via a reformulation of ( |22] | in terms of 
ADM variables on the slice iQ, 



where 



i^jk- 



(24) 



(25) 



The identification of the Newman-Penrose ^'4 with the grav- 
itational radiation content of the spacetime is a result of the 
peeling theorem, which states that in an appropriate frame the 
^'4 component of the curvature has the slowest falloff with ra- 
dius, 0(l/r). The conditions of this theorem are not satisfied 
exactly at a small radius and in the chosen frame. While there 
are proposals for how this situation can be improved ll44ll . we 
find that beyond r-^ > 30 M in fact our measure of '^4 scales 
extremely well with the different extraction radii r-^ , suggest- 
ing that the peeling property is satisfied to a reasonable ap- 
proximation (see Fig. nil. 

The gravitational wave polarization amplitudes /i+ and hx 
are related to ^'4 by i45J 



(26) 



where the double over-dot stands for second-order time 
derivative. The flux of linear momentum emitted in gravi- 
tational waves in the i-direction can be computed from the 
Isaacson's energy-momentum tensor and can be written in 
terms of the two polarization amplitudes as Jstl 



IGtt 



dfl 



hi 



(27) 



where rii = Xi/r is the unit radial vector that points from the 
source to the observer and dil = sin 9d(pd9 is the line element 
of our extraction 2-sphere S^. Using Eq. ( |26] |. this leads to 
an expression for the momentum flux in terms of ^4 as it is 
commonly used in recent numerical relativity calculations |0, 

filHIiililliilii: 



J-j — lim 



'sch 

167r 



dn ' 



dt*4 



(28) 



The Schwarzschild radius, rgch, is derived from the coordinate 
(isotropic) radius via the standard formula 



^sch 



1 - 



M 



2ru 



(29) 



assuming a constant ADM mass M = ^^adm throughout 
the simulation. With this choice of radial coordinate, expres- 
sion ( |28] | has been shown to provide recoil velocities which 
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are in better agreement with those obtained through gauge- 
invariant perturbations than with the ahernative coordinate ra- 
dius, {cf. Sect. IIII Bb and reported in the literature (Additional 
details on the numerical measurement of ^'4 are presented in 
AppendixlBl) 

B. Kick measurements via gauge-invariant perturbations 

An independent method to compute the linear momentum 
carried away by gravitational radiation is based on the mea- 
surements of the non-spherical gauge-invariarit perturbations 
of a Schwarzschild black hole (see Refs. ifsol Isii I52I1 for ap- 
plications to Cartesian coordinates grids). In practice, a set of 
"observers" is placed on 2-spheres of fixed coordinate radius 
Tj, , where they extract the gauge-invariant, odd-parity (or ax- 
ial) current multipoles and even-parity (or polar) mass 
multipoles Q^^^ of the metric perturbation ifssll . The numer- 
ical implementations of the gauge-invariant variables is done 
by following the multipolar analysis outlined by Abrahams 
and Price ||5J|. The and Q^^^ variables are related to /i+ 
and hx as llssll 




Here _2Y^'^ are the s = —2 spin-weighted spherical harmon- 
ics and {£, m) are the indices of the angular decomposition. 
Validations of this approach in 3D vacuum spacetimes can be 
found in Refs. lf52ll5a,l57[| . while its use with matter sources 
has first been reported in 115811 . 

We note that the notation introduced in Eq. dSOl l could be 
misleading as it seems to suggest that is, always of odd- 
parity and hj^ is always of even-parity. Indeed this is not 
true in general and in the absence of axisymmetry, i.e., when 
TO 7^ 0, both hx and h+ are a superposition of odd and 
even parity modes. It is only for axisymmetric systems, for 



which only to = modes are present, that Q^,-^ and Q^^^ 
are real numbers, that h+ is only even-parity and /ix is only 
odd-parity. Despite this possible confusion, we here prefer to 
maintain the notation of Eq. dSOb which is the most common 
in the literature ifssll . 

The flux of linear momentum emitted in gravitational waves 
in terms of Q^^^ and can be computed by inserting 
Eq. ( [30] l in Eq. ( |27] ). then decomposing rii in spherical har- 
monics and performing the angular integral. This proce- 
dure goes along the lines discussed by Thorne in Ref. lf59ll . 
where all the relevant formulae are essentially available {cf. 
Eq. (4.20) there. See also Ref. ll60ll 1. so that we only need to 
adapt them to our notation. In Ref. lf59tl the even-parity (or 
electric) multipoles are indicated with and the odd-parity 
(or magnetic) ones with Sim- They are related to our notation 
by 

= , (31) 
= Q"^^ , (32) 

where '^^\fim = d^fem/dt^. From the well known property 
(Qtm^)* ~ (~-'^)™Qf'-m' where the asterisk indicates com- 
plex conjugation, one can rewrite Eq. (4.20) of Ref. (5^ in 
a more compact form. Following Ref. 1I6II1 where the lowest 
multipolar contribution was explicitly computed in this way, 
it is convenient to combine the components of the linear mo- 
mentum flux in the equatorial plane in a complex number as 
+ iJ^y The multipolar expansion of the flux vector can be 
written as 

00 t 

Tx+iTy^Y.Y. (•^''" + i-^D ' (33) 

1=2 m=0 
00 t 

^z = T.T. ' (34) 

i=2 m=0 

where Sm = 1 if to 7^ and (5,„ = 1/2 if to = 0. Each 
multipole reads 



2i 



+ n+ nx 



-niQem-l + '^ImQtmQi -(m+1) 



y 167r^(£+l) 

K (0+ Q+ Qf^, ,)+b+ (Qt Q+ , v^.^, , 

^™ — m^£+l m— 1 ' ^£ — + 1 m— 1 j ' tm \ ^£m^£+l —(m+1) ^£m^£+l —(m+1) 



(2^+ l)(2£ + 3) _ 



— ( '^y 



?7r£(£+ 1) 



2to Im 



(2^+ l)(2£ + 3) 



Re 



1 1 



(35) 
(36) 



and Note that here both jFf;™ and JF,, are real numbers and are 



a 



± 



^/{^±m){(.Tm+l) , (37) 
bf„^ = v/(^±TO+l)(£±m + 2) , (38) 
ctrn = x/ie-m + l){£~m+l) . (39) 
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obtained as the real and imaginary part of the right-hand-side 
of Eq. ( [35] l. For a general system without symmetries one is 
expecting jFf™ to be nonzero. However, our initial data set- 
up, an inspiraling binary with spins anti-aligned and parallel 
to the orbital angular momentum, implies that the linear mo- 
mentum flux vector is completely contained in the equatorial 
plane of the system and so that JF|"' = by construction. 
Since we are imposing equatorial symmetry {i.e., invariance 
for IT 6) we have that multipoles with £ + m ^ even 
are purely even-parity (i.e., Qf^^ ^ and Q^^^ = ) and 
those with ^ + m = odd are purely odd-parity (i.e., Q\„^ = 
and Q^^j ^ 0). As a final remark, we note that for £ = m ^ 2, 
our Eq. ^ reduces to Eq. (9) of Ref. lHH]. 

IV. RESULTS 

This section collects the results of our analysis of the re- 
coil velocity of spin-aligned binaries and discusses the differ- 
ent aspects of the study which combined provide a consistent 
and accurate picture of this process. We will first concentrate 
on the systematic error introduced by the use of initial data 
with zero Unear momentum and on the techniques we have 
developed to remove it. We will then discuss the actual com- 
putation of the recoil velocities and their dependence on the 
spin ratio, highlighting the modes of the radiation which are 
largely responsible for the asymmetric emission. Finally, we 
will discuss the accuracy of our measurements and our ability 
to preserve mass and angular momentum to below 1%. 




t (M) 



FIG. 5: The recoil velocity of the binary r O is compared to those of 
the same system but with either a larger or a smaller initial separation 
{i.e., rOl and rOs, respectively). Note the same recoil velocity is 
obtained when the integration constant is properly taken into account, 
while an error as large as ~ 13% is made otherwise. 



A. Initial transients in the waveforms 

Both Eqs. (l28T l and ( l35T l provide an expression of the re- 
coil velocity in terms of the radiated (linear) momentum per 
(infinitesimal) time interval. A time-integration of those equa- 
tions is needed in order to compute the recoil and this obvi- 
ously opens the question of determining an integration con- 
stant which is in practice a vector. Fortunately, this integration 
constant has here a clear physical meaning and it is therefore 
easy to compute. In essence it reflects the fact that at the time 
the simulation is started, the binary system has already accu- 
mulated a non-vanishing net momentum as a result of the slow 
inspiral from an infinite separation. 

Since the initial data is constructed so as to have a vanishing 
linear momentum, there will be a inconsistency between this 
assumption and the actual evolution of the initial data. Stated 
differently, the numerical evolution of the Einstein equations 
will soon tend to a spacetime which is different from the ini- 
tial one and indeed corresponding to one with a net linear 
momentum. This momentum is the one that the binary has 
gained when inspiralling from t = —oo till t = 0. Calculat- 
ing the integration constant amounts therefore to computing 
the vector accounting for this mismatch and is essential for 
a correct measurement of the recoil velocity. The error made 
when neglecting this constant, as routinely done in numerical- 
relativity calculations, inevitably produces a systematic devi- 
ation from the correct answer and, as we will show in the next 



section, it can altogether prevent from having even the quali- 
tative behavior right. 

The relevance of this integration constant depends on the 
initial separation and it is more important for binaries that 
start their evolution akeady quite close. This is rather obvi- 
ous: the tighter the binary is, the larger the emitted momen- 
tum per unit time and the more important is to evaluate the 
initial mismatch. Fig. |5] helps to illustrate this point and can 
be discussed before entering into the details of how we actu- 
ally compute the integration constant. The figure shows the 

time evolution of the recoil velocity |w|kick = y^f ^ + for 

the same binary system having spin ratio ai/a2 = — 1 but 
with increasing initial separation. More precisely, we consider 
systems rOl, rO and rOs which differ only in the initial sepa- 
ration, which is about 8.4, 6.0 and 5.6 M, respectively. The 
data Fig. |5]is properly shifted in time so as to have the curves 
overlap and shows that only when the integration constant is 
properly taken into account, do the three simulations yield the 
same recoil velocity [cf., solid, dashed, and dotted lines). On 
the other hand, when the integration constant is not included 
in the calculation, different evolutions will yield different esti- 
mates, with a systematic error that can be as large as 13% (cf., 
long-dashed and dot-dashed lines) and is clearly unacceptable 
given that the overall precision of the simulations is below 1% 
[cf.. Figs. [TTI - [T2l and the discussion in Sect. lIV Dt . 

Besides providing the right answer, the calculation of the 
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FIG. 6: Left panel: Evolution in velocity space of the recoil-velocity vector. Very little variation is recorded before the radiation reaches the 
observer at r^, = 50 M (dotted lines in the two insets). The absence of the proper linear momentum in the initial data triggers a rapid and an 
almost straight-line motion (dashed line) of the center of the spiral away from the origin of coordinates during the initial stages of the evolution. 
After this transient motion, the evolution is slower, with the spiral progressively opening up (solid line). The vector to the center of the spiral 
corresponds to the initial linear momentum of the spacetime and is used as integration constant for Eqs. l l28b and ( I35t . The final part of the 
evolution is characterized by a change in the spiral pattern (long-dashed line) as a result of the interaction of different modes in the ringdown 
of the final black hole. Note that the figure has been rotated clockwise of about 30° to allow for the two insets. Right panel: Initial behavior of 
the recoil velocity (upper graph) and of the waveform {Q22) for model rO (lower graph). This figure should be compared with the initial vector 
evolution of the recoil velocity shown in the left panel where the same types of lines have been used for the different stages of the evolution. 
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FIG. 7: Left panel: The same as in the left panel of Fig.|6]but for system rl. Shown in the inset is the sudden re-orientation of the recoil 
velocity vector during ringdown and corresponding to a new spiral with different aperture (long-dashed line). Although more pronounced in 
r7, the appearance of this "hook" at ringdown is seen all the members of the sequence. Right panel: The same as in the left panel of Fig.|6] 
but for system r7. The upper graph concentrates on the final stages of the evolution in of the recoil velocity and on the appearance of a second 
peak during ringdown (long-dashed curve). The lower graph shows the same but in terms of the Q22 waveform. A discussion of these final 
stages of the evolution is made in Sect. llV Cl 
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integration constant also results in a considerable saving in 
computational costs. The complete dynamics of the binary 
rOl including the merger and ringdown, in fact, requires sim- 
ulations for about 600 AI; the same answer in terms of recoil 
velocity can be obtained with the system rOs, whose dynam- 
ics is fully accounted for with a simulation lasting only for 
340 M. 

Having stressed the importance of including the integra- 
tion constant in the measurement of the recoil velocity, we 
next illustrate how to actually compute it. In essence, it is 
sufficient to look carefully at the evolution in the velocity- 
space of the two components Vx and Vy of the recoil veloc- 
ity (because of the symmetry the z-component is zero but the 
method described here can be easily extended to the case in 
which 0). This is shown in the left panel of Fig. |6l 

which reports the track of the "center of mass" for system 
rO in such a space. Different types of line refer to differ- 
ent intervals in time during the evolution and, for an ob- 
server at — 50 Af , the dotted one refers to t < 50 M, 
the dashed one to 50 M < i < 75 M, the continuous one 
to 75 M < t < 183 M, and finally the long-dashed one to 
t > 183 M. 

Clearly, for t < 50 M the system undergoes very little evo- 
lution in velocity-space (c/, dotted line in the inset within the 
inset of the left panel) but a rapid change, lasting for about 
25 M, takes place as the radiation reaches the observer The 
radiation received has information about the "correct" linear 
momentum of the spacetime which is solution of the Einstein 
equations for system 7'0 as if it had inspiraled from infinity, 
and thus rapidly moves the center of mass to a net nonzero 
recoil velocity (cf., almost-straight dashed line in the inset in 
the left panel). Once the system has adjusted for the proper 
linear momentum, the evolution proceeds as expected, with 
the recoil velocity vector slowly tracking a spiral in velocity 
space. This is an important point which we prefer to under- 
line: the rate of change of linear momentum is very large only 
initially and this is because as the binary migrates from the 
initial non-radiating state (the data is conformally flat) to the 
consistent radiating state, it will emit the amount of linear mo- 
mentum it would have emitted when inspiralling from infinite 
separation. After this burst of linear momentum, the evolution 
of the recoil velocity is minute, essentially until it grows very 
rapidly during the last orbit. 

Computing the integration constant consists then in cal- 
culating the position of the center of the spiral and this can 
be done either by a simple inspection of a graph in the 
velocity-space, from which compute the center of the spiral 
or, equivalently, by searching for the initial vector that would 
lead to an essentially monotonic in time growth of the recoil 
velocity It^]. The latter procedure does not require a human 
judgment but we have found it to yield the same answer (to 
less than 1 km/s) as the one guessed by looking at the veloc- 
ity space. 

The right panel of Fig. |6] shows the same evolution as the 
left one, but through different quantities. The upper panel, 
in particular, shows the time evolution of the recoil velocity 
and the rapid changes it undergoes initially when the radia- 
tion first invests the observer. The lower panel, on the other 



hand, shows the Q22 amplitude and highlights that, while the 
initial burst of radiation stops after t ~ 50 M {cf., dotted line), 
the waveform is still not fully consistent until t ^ 75 AI (cf., 
dashed line). 

The procedure discussed so far for the calculation of the 
integration constant relative to the binary system rO applies 
qualitatively to all the other members of the sequence, with 
differences that are due essentially to the times at which the 
various stages take place. 

It is worth remarking that the evolution of the recoil vector 
in the velocity-space has another interesting feature during the 
final stages of the evolution and when the final black hole is 
ringing down. This is marked as a long-dashed line in the left 
panel of of Fig.|6]and shows a break in the building of the spi- 
ral and the appearance of a new spiral with a different aperture 
(we refer to this feature as "the hook"). This is more evident 
in the left panel of Fig.|7] which shows the evolution of the re- 
coil vector for the binary system r7 and offers a magnification 
of the hook in the inset. A more detailed description of this 
feature is beyond the scope of this paper and will be presented 
in a future work, but we can here point out that the hook ac- 
counts for a rapid change in the recoil velocity and it is due to 
the interplay of different modes during the ringdown. This is 
clearly illustrated in the right panel of Fig. |7] which similarly 
reports the time evolution of the recoil velocity and the final 
stages of the Q22 waveform. 



B. Recoil velocities 

The recoil velocity has been calculated for the sequence of 
models listed in Table U As mentioned in Sect. IllDI this se- 
quence corresponds to equal-mass black holes, whose initial 
spins are unequal, though always aligned with the z-axis. The 
rO model has equal but opposite spins, while the r8 model 
has equal and aligned spins on the black holes, with other 
models corresponding to intermediate values, as outlined in 
Section HTD] Since the total initial orbital angular momentum 
L of the system is chosen to be constant over the sequence, 
the initial separations of the black holes increases in the se- 
quence, as well as the time to merger due to spin-spin effects 
which contribute to an orbital "hang-up" in the aligned case. 

We extract gravitational waves by both the gauge-invariant 
and the 'I'4 methods described in the previous section and by 
interpolating the radiation-related quantities onto 2-spheres at 
coordinate radii = 30 M, 40 M, 50 M, and 60 AI. The use 
of multiple extraction radii is made to check the consistency 
of the measurement and the precise value of the extraction 
radius has little influence on the actual kick calculation. In the 
case of the binary system tQ we have verified that the recoil 
velocity yields the same value with differences that are smaller 
than 2 km/s for extraction 2-spheres at distances larger than 
30 M. As a result, we have used = 50 AI as the fiducial 
distance for an observer in the wave-zone and all of the results 
presented hereafter will be made at this extraction 2-sphere. A 
validation that the gauge-invariant quantities have the proper 
scaling with radius is presented in AppendixICl 

The evolution of the recoil velocity for the entire sequence 
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FIG. 8: Left panel: Recoil velocity as a function of the spin asymmetry parameter axja^ for the models listed in Table |l] Indicated with a 
continuous lines are the results obtained via ^'4, while a dashed line is used for the gauge-invariant quantities ■ Right panel: Final recoil 
velocity calculated with both the use ^4 (empty circles) and the gauge-invariant quantities (stars). Shown in the inset is the incorrect scaling 
obtained when the correction for the integration constant is not made. 



listed in Table U is displayed in the left panel of Fig. |8] It is 
apparent that the suitable choice of the integration constant 
discussed in the previous section yields early evolutions that 
are always monotonic in time and that, as expected, the largest 
recoil velocity is generated for the case in which the asymme- 
try is the largest, namely for the binary rO. The left panel 
Fig. |8] also shows that the profile for each case is rather sim- 
ilar, with the largest contribution to the kick velocity being 
generated in a period of about 80 M, corresponding roughly 
to the timescale of the last orbit and merger. Furthermore, it is 
notable that 95% of the acceleration occurs ~ 30 M after the 
appearance of the first common apparent horizon, indicating 
that the kick is generated not only by the final stages of the in- 
spiral (i.e., by the "plunge") but also and more significantly by 
the ringdown of the final black hole. This fact helps to explain 
why accurate recoil velocities can be obtained by evolutions 
involving very few cycles only, provided the integration con- 
stant is properly taken into account. 

It is worth noting that during the final stages of the evo- 
lution, the recoil velocity is not monotonic but shows at least 
two peaks, whose relative amplitude depends on the spin ratio. 
For spin ratios ^ — 1 the first peak is hardly visible, while the 
second one is the most pronounced one. As the spin ratio in- 
creases, however, the first peak becomes more prominent and 
for spin ratios ~ 1 it becomes comparable with the second 
one or even larger for binaries rQ and r7. As mentioned in the 
previous Section and further discussed in the following one, 
the appearance of these peaks is related to the interplay of dif- 
ferent mode-contributions during the ringdown. The second 



peak, in particular, can be associated to a rapid change in the 
recoil-velocity vector and is behind the characteristic "hook" 
discussed in the left panels of Figs. |6] and [T] While additional 
work is needed, especially in thorough perturbative investi- 
gations, to fully account for the rich, post-merger properties 
of the recoil velocities, we believe the double-peak evolution 
to be physically genuine since it is seen in all binaries and is 
supported by the highly accurate and convergent simulations. 
As a representative measure of the accuracy in determining 
these recoil velocities, we mention that we have carried out 
simulations also for the binary system r8, in which the black 
holes have identical spin and thus from which no kick should 
result. The computed recoil velocity has been found to be 
10~^ km/s, clearly indicating that our evolutions do an excel- 
lent job in preserving the orbital symmetry of these binaries. 

We have found that the evolution of the recoil velocity gen- 
erated by spin asymmetries appears to be rather different from 
the one generated by mass asymmetries It8|, lifi, i47il and which 
shows much larger variations between the maximum attained 
value and the final one. Once again, this different behavior is 
related to the different interplay of the ringdown modes in the 
case of mass asymmetries and will be presented in a separate 
work. 

The recoil velocities attained by the final black holes and 
shown for in the left panel of Fig |8] can be studied in terms of 
their dependence on the spin ratio 01/02, which can also be 
regarded as the "asymmetry" parameter of the system, being 
the largest for 01/02 = —1 and zero for 01/02 = 1. These 
velocities are collected in Tablelllland are shown as a function 



14 



TABLE 11: Final kick velocities in units of km /s for the models listed 
in Tab. |I] Columns two and three show the values obtained using 
the gauge-invariant quantities Q^^]^ and 'I'4 respectively and taking 
into account the integration constant. Columns four and five, on the 
other hand, show the results obtained when ignoring the integration 
constant. The same data are shown in the right panel of Fig. [8] 



Model 




*4 




\l/4, no ic 


rO 


263.2 


261.8 


288.9 


288.4 


rl 


222.4 


221.4 


211.9 


210.6 


r2 


187.1 


186.2 


174.8 


173.3 


r3 


143.3 


144.0 


155.9 


157.3 


r4 


104.8 


106.1 


100.0 


101.3 


r5 


81.4 


81.5 


76.9 


77.0 


r6 


45.6 


45.9 


55.4 


56.2 


r7 


19.4 


20.6 


13.8 


14.8 


r8 


0.0 


0.0 


0.0 


0.0 



of ai/a2 in the right panel of Fig [8] where we have indicated 
with open circles the values obtained using ^'4 and with stars 
those obtained using the gauge-invariant perturbations. 

The data in the right panel of Fig|8]is shown together with 
its error-bars, which include errors from the determination of 
the integration constants, from the truncation error and from 
the amount of ellipticity contained in the initial data. We have 
estimated these errors to be of 5 km/s for binaries r0-r5 and 
of 8 km/s for binaries r6 and rl. Shown also in the inset is the 
recoil data obtained when ignoring the integration constant. It 
is remarkable that when the proper evaluation of the initial 
transient is not made, the data does not show the remarkable 
correlation with the spin ratio which is instead shown by the 
corrected data. Quite surprisingly, however, the correlation 
found the one predicted by PN studies. We recall, in fact, that 
using PN theory at the 2.5 order, Kidder ll62ll has concluded 
that in the case of a circular, non-precessing orbit, the total 
kick for a binary system of arbitrary mass and spin ratio can 
be expressed as jst] 



l^'lkick = cr 



(1 



C2- 



= £202 1 - 



Ol 
0,2 



(40) 



where q = M1/M2 is the mass ratio and is equal to one for 
the binaries considered here, thus leading to the second form 
of Eq. ( l40b . The coefficients ci and 02 = 02/ 32 depend on the 
total mass of the system and on the orbital separation at which 
the system stops radiating, which is intrinsically difficult to 
determine with precision since it lies in a region where the PN 
approximation is not very accurate. Indeed, we find that the 
coefficient C2 is not really a constant in the case of equal-mass 
binaries but, rather, it can be seen to depend at least linearly 
on the spin ratio. 

This is shown in Fig.lH whose upper panel offers a compar- 
ison among the computed data for the recoil velocity (open 
circles) with the least-squares fits using either a linear (dot- 
ted line) or a quadratic dependence (dashed line). It is quite 
apparent that a linear dependence on 01/02, such as the one 



250 
200 



I 150 



.2 100 



50 



-20 - 



-e- 



data 
Im. fit 
quad. fit_ 




FIG. 9: Upper panel: Comparison of the computed data for the re- 
coil velocity (open circles) with the least-squares fits using either a 
linear (dotted line) or a quadratic dependence (dashed line). Lower 
panel: Point-wise residuals computed with the linear (dotted line) or 
a quadratic fit (dashed line). 



expected in Eq. (|40] | for C2 = const, does not reproduce well 
the numerical data and yields point-wise residuals of the order 
of 20 km/s. These are shown with a dotted line in the lower 
panel of Fig. |9] A quadratic dependence on ai/02, on the 
other hand, reproduces the numerical data very nicely, with 
residuals that are of the order of 5 km/s, as shown with a 
dashed line in the lower panel of the same figure, and thus 
compatible with the reported error-bars. 

We can re-express Eq. ( |40l i in the more generic form 



|u|kick 02 



Ol 
02 



l«2|/ 



(41) 



where 02 plays here the role of a "scale-factor". The function 
/(01/02) with 01/02 e [—1, 1] and maximum at 01/02 = 
— 1 can then be seen as to be determined from numerical- 
relativity calculations (or higher-order PN approximations) 
and our least-squares fit suggests the expression 



/quad. - 109.3-132.5 



23.1 



km/s . 
(42) 



The maximum kick velocity for a given 02 is then readily 
calculated even without a detailed knowledge of the function 

/(01/02) as 



(k|kickr^^(02) = |02|/(-l). 



(43) 
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FIG. 10: The total kick calculated via Eq. ( 135b up to ^ = 7 is compared to the contributions of individual terms qi and §2, as well as the sum 
of term excluding these. In the case of the rO system (left panel) the spins are anti-aligned and the g2 term is dominant and the qi term does 
not provide a significant contribution. In the case of the r7 system (right panel), on the other hand, the spins are essentially aligned and the 
while the q2 term is still dominant, the qi term also makes a significant contribution. 



Using the data reported in Tablellllfor a-? = —0.584 we obtain 
for 1 02 1 = 1 that the maximum recoil velocity attainable from 
a binary system of equal-mass black holes with spins aligned 
to the orbital angular momentum is 448 ± 5km/s. This is 
in very good agreement with our previous estimate made in 
Ref. lUOtl with a smaller sequence and in equally good agree- 
ment with the results reported in Ref. i^. 

C. Mode contributions to the recoil velocity 

For the models studied in the previous section we have eval- 
uated Eq. ( |35] ) including modes up to ^ = 7. In practice, 
however, we find that the recoil is strongly determined by the 
lower-mode contributions. In particular, the two terms 

1 /"so" • 

'Zi^^VyQ^sQa^-s, (44) 
q^ = -^Q+_^Q-^ (45) 

are the dominant ones. This can be seen in Fig.[TO] where the 
time evolutions of the terms qi and q2 are plotted (dotted and 
dashed lines, respectively) together with the total kick calcu- 
lated via Eq. (l35T l (solid line), and with the contributions from 
all other terms up to £ = 7 excluding qi and q2 (long-dashed 
line). A rapid inspection of the figure reveals that the kick is 
dominated in particular by the q2 term, whereas the qi term 



has a magnitude of the order of all the other modes combined. 
A similar result holds for each member of the sequence, so 
that the two contributions determine the final kick to more 
that 95%. It should be noted that the mode contributions are 
vector quantities, just as the kick velocity itself, and are not 
always aligned or even maintain the same angle to each other 
during the duration of the recoil. 

This coupling also goes some way to explain some features 
of the recoil velocity profiles displayed in Fig. [8] As men- 
tioned in the previous section, in fact, the binaries r4 to r8 
show a clear double peak in the evolution of the kick velocity 
before it settles down to the final value. The same feature can 
also be seen in the more asymmetric rO to r3 binaries, where 
it appears as a flattening of the slope near the maximum. Since 
the two peaks are shown both by the gauge-invariant and by 
the ^'4-based techniques (which are rather different in both 
the assumptions they rely on and in the practical implementa- 
tion) we do not believe them to be a simple numerical artifact. 
Overall, the properties of the recoil velocity near its maxi- 
mum, and before it settles to the final value, are determined by 
the relative phases of the two contributions identified above. 
An analysis of the terms qi and q2 in vector-space, and which 
will be presented in a subsequent work, reveals that when they 
are relatively aligned at the peak of the acceleration, there is 
a clear single peak in the evolution. For the more symmet- 
ric models, on the other hand, the two contributions are more 
anti-aligned and a double peak results. 

These considerations in the vector evolution of the two con- 
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tributions qi and q2 need also to be linked with the evolution in 
vector space of the recoil velocity. As stressed in Sect. |1V A| 
in fact, there is a distinct kink in the evolution of the veloc- 
ity vector towards the final stages of the merger (this feature 
is indicated with a long-dashed line in the Vx vs. Vy plots of 
Figs. |6] and I?]). The presence of the kink corresponds to a lo- 
cal decrease of the recoil velocity and hence to the minimum 
between the two peaks. Because this decrease is more pro- 
nounced for the lower-kick binaries r4 to ?'8, the first peak 
becomes more evident there. 



The second method instead, assumes that the final black 
hole has settled to a Kerr one and uses the the rotational- 
induced distortion of the apparent horizon of the final black 
hole to estimate its spin. Defining Cp and Cg to be respec- 
tively the apparent horizon's polar and equatorial proper cir- 
cumferences, their ratio Cr = Cp/Ce will undergo damped 
oscillations as the perturbed black hole settles to a Kerr state 
through the quasi-normal ringing. The final value of Cr can 
be expressed as a nonlinear function of the dimensionless spin 



parameter a J/M^ as ll2Ul66ll67ll6: 



D. Angular Momentum and Mass Conservation 

In this section we discuss the radiated angular momentum 
and energy during the evolution of the different initial-data 
sets. We compute the radiated angular momentum and mass 
by calculating their difference between the initial data and that 
of the final black hole, and then compare these quantities with 
the corresponding ones measured in terms of the emitted grav- 
itational radiation. The differences in the two independent es- 
timates serve therefore as stringent indicators of the conserva- 
tive properties of our code. 

The radiated angular momentum can be simply written as 
the difference between the initial and final values 



' rad 



J fin Jh 



(46) 



where, as a result of the conformal flatness of the initial-data 
slice, Jini is given by the simple expression (see for exam- 
ple |l39l|M.|41D and discussion in Sect.|llll 



J. 



Ci X + C2 X + Si + S2 ■ (47) 



Here Ci, p^ and Si are the position, the linear momentum and 
the spin of the i-th black hole. The final angular momentum 
Jfin, on the other hand, is set to be equal to the spin of the final 
black hole after all the radiation has left the computational do- 
main. Two different methods are used to obtain this measure, 
both of which are based on properties of the apparent horizon 
of the final hole. 

The first method employs the isolated/dynamical horizon 
formalism and searches for a rotational Killing vector on 
the final apparent horizon so as to measure the spin of the final 
black hole as 1631 



Kab 



(48) 



We note that this expression ( |48] l is valid on any sphere where 
a Killing vector 0" can be found, and is therefore a quasi-local 
measure of the angular momentum. In particular, at large dis- 
tances where the spacetime is close to axisymmetric, there is 
a good approximation to an angular Killing vector, and we 
can apply this expression to determine the angular momen- 
tum of the spacetime. Note also that Eq. ( |48] | is identical to 
the ADM angular momentum when evaluated at spacelike in- 
finity. (Refs. ||64[ I65I1 also give a quasi-local formula for the 
angular momentum flux due to gravitational radiation.) 



a (a) 



1 + vT" 



-E - 



(49) 



where E{k) is the complete elliptic integral of the second kind 



.7^/2 

m = / 

Jo 



Vi - fcsin^ ede. 



(50) 



By inverting numerically Eq. (I49t we obtain a from the late 
time Cr that is measured from the apparent horizon shape. 
Note that for computing J we need to multiply a by the square 
of the final mass, which we take to be M^j^^j — Miad- An al- 
ternative choice involving the total mass Eq. (fTTI i as measured 
from the apparent horizon would lead to essentially the same 
results. 

As mentioned at the beginning of this section, the determi- 
nation of the radiated angular momentum can also be done 
using directly the asymptotic waveform amplitudes /i+ and 
/ix as lUliilTSl 



dtdh 



(51) 



where the amplitude /i+ and themselves can be expressed 
either in terms of the Zerilli-Moncrief gauge-invariant vari- 
ables (5^„, Q^^rn alternatively, in terms of the Newman- 
Penrose scalar *I'4. A comparison between the two approaches 
is presented in Appendix |Cl where it is shown that the dif- 
ferences are minute. Because of this, hereafter we will re- 
fer to asymptotic amplitudes measured in terms of the gauge- 
invariant variables only. Additional details on the resolution 
of the extraction 2-sphere are also presented in AppendixlB] 

The left panel of Fig. [TT] summarizes this comparison by 
showing, as functions of the spin ratio ai/a2, Jun from 
Eq. ( |48] |. Jiad from Eq. ( BTl i both adding nicely to yield 
Jini. Note that Jini is growing linearly as it is obvious from 
Eq. ( |47] ), but also that that a similar behavior is shown by the 
radiated angular momentum (and hence by the final spin of the 
black hole). Using a linear fitting we can derive phenomeno- 
logical expressions for the relative losses of angular momen- 
tum 



~T '^rad I 



Xrad 



and the relative spin-up of the final black hole 



■/fin _ ^,7 / Oi 

— — Sfin I 

Jini V ^2 



Xfin 



(52) 



(53) 
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FIG. 11: Left panel: Dependence on the spin ratio of the initial total angular momentum Jini [as computed from Eq. l|47H, of the radiated 
angular momentum Jrad [as computed through the gauge-invariant waveforms], and of the final spin of the black hole Jfln- All quantities 
show a linear behavior, whose coefficient are collected in Table HIH Right panel: Relative error A J/ Jini in the conservation of the angular 
momentum [cf., Eq. ||54H. Different curves refer to whether the final spin of the black hole is computed using the isolated/dynamical horizon 
formalism (triangles) or the distortion of the apparent horizon (squares). In both cases the error is of about 1% at most for simulations at the 
medium resolution. 




FIG. 12: Left panel: Dependence on the spin ratio of the ADM mass Mj^^^-^, of the scaled radiated energy M^ad [as computed through the 
gauge-invariant waveforms and scaled by a factor of 10 to make it visible], and of the final mass of the black hole Affln. All quantities show 
linear behaviors, whose coefficients are collected in Table HiH Right panel: Relative error Ai\//Mini in the conservation of the energy \cf., 
Eq. \5(M. Note that the error is of about 0.5% at most for simulations at the medium resolution. 
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TABLE III: Coefficients for the phenomenological expressions l l52t 
and l l53t (and tlie corresponding coefficients for AA/rad, fln/A^) by 
means of wfiich it is possible to compute the relative losses of energy 
and angular momentum, as well as the final mass and spin of the 
black hole in binary mergers in which the spins are orthogonal to the 
orbital plane. 



A-' 

Srad 
Xrad 


0.0513 
0.2967 


Srad 
Xrad 


0.0118 
0.0437 


Xfln 


-0.0513 
0.7033 


M 
Xfln 


-0.0118 
0.9563 



The fitted values for ^;4d fin ™^ x/ad fin presented in Ta- 
ble Hn] and readily indicate that the system looses 24% of its 
initial orbital angular momentum in the case of anti-aligned 
spins and up to 34% for aligned spins. 

To the best of our knowledge expressions ( |52] i and ( |53] l do 
not have a PN counterpart and yet, since they depend only 
on the spin-ratio, they represent simple and powerful ways 
of estimating both the efficiency in the extraction of angular 
momentum and the spin of the final black in a binary merger 
when the spins are orthogonal to the orbital plane. This infor- 
mation could be easily injected in those A^-body simulations 
in which the interaction of binary black holes is taken into 
account iFtiIi and thus yield accurate estimates on final distri- 
bution of black-hole spins. 

Since we have two independent and different ways of com- 
puting Jiad [i.e., either from Eq. ( fSTT i or from Eq. ( |46] ll we can 
quantify our ability to conserve angular momentum by mea- 
suring the normalized residual 

^fin ^" ^rad ^ini (^A\ 

'-'ini 'Jini 

This is shown in the right panel of Fig. [TT|and the two differ- 
ent lines refer to the two measures of the final spin of the black 
hole, i.e., either via the isolated-horizon formalism (triangles) 
or via the distortion of the apparent horizon (squares). In both 
cases the error is extremely small, ranging between 1.1% and 
0.2% for simulations at the medium resolution, and thus pro- 
viding convincing evidence of our accuracy in the preserva- 
tion of angular momentum. It should be noted that while there 
seems to be a small advantage in using the isolated horizon 
measure, the differences are too small to be significant. In- 
deed, a small change in the procedure, such as the use of the 
mass measured via the apparent horizon via Eq. ( |49l ) in place 
of Mini — Mfin (as we are doing in this figure), would revert 
the advantage. 

We proceed next to a similar analysis for the conservation 
of the mass-energy of the system by considering the difference 
between the the initial mass and final plus the radiated masses. 
As for the initial mass we obviously consider the ADM mass 
of the system Af^^^j, while the radiated energy Mrad is com- 
puted through the gravitational waveforms lf55ll72l PTSll 



dtdVt 



r- 
16^ 



(55) 



functions and to use as final mass of the black hole ilffi„, the 
one given by Eq. ( fTTl ) and measured via the apparent horizon. 

The left panel of Fig. [T2l shows A/^^^,, M^n and Afrad, 
with the latter rescaled the radiated by a factor of ten to make 
it more visible. Also in this case there is a clear linear behavior 
of both the radiated energy and of the final mass of the black 
hole in terms of the spin ratio. As a result, phenomenologi- 
cal expressions of the type ( l52b and ( l53T l are possible also for 
A/fi„ and A/iad- The corresponding values of the coefficients 
Crad, fin Xrad. fin also presented in Table Hill 

Finally, to check the precision at which the energy is con- 
served, and in analogy to Eq. ( l54b . we have computed the rel- 
ative error 



AAf ^ A/fin + Afrad - Mj, 



(56) 



As for the angular momenta, we have chosen to express the 
right hand side of Eq. dSSl l in terms of the Zerilli-Moncrief 



and plotted this as a function of the spin ratio in the right panel 
of Fig. [T2I Clearly, also the energy losses are extremely small 
and for all the binaries in the sequence, the error in the energy 
balance is below 0.52% at the medium resolution. Table HVl 
summarizes the numerical results for the radiated energy and 
angular momentum for the members of the sequence. 



V. CONCLUSIONS 

We have performed a highly-accurate study of recoil veloc- 
ities in binary black hole mergers from a sequence of equal- 
mass black holes with varying spin configurations. In this se- 
quence, the spins are aligned with the orbital angular momen- 
tum since there are strong indications that such alignment is 
preferred in astrophysical situations. This makes our choice 
of initial data especially realistic and our results particularly 
relevant also within an astrophysical context. 

In practice, the initial configurations are built so that the 
spin of one of the black holes is kept at a constant dimension- 
less value a2 = 0.584 while the other varies from ai = —02 
to fli = +a2, thus spanning a range between —1 and 1 in spin 
ratio. We have followed our black hole evolutions for about 
two to four orbits and then throughout the plunge, merger, 
and ringdown phases. This work thus extends and refines re- 
cent results obtained from a reduced but similar initial-data 
sequence IToll . 

The main aspects of this work, which revolve around the 
methods used, the tests performed and the results obtained, 
can be summarized as follows. 

Methods. To increase the significance of our results and our 
confidence in their accuracy, we have implemented two inde- 
pendent methods for the calculation of the linear momentum 
from the emitted gravitational radiation. These are based on 
either the measure the Newman-Penrose scalar or on the 
calculation of the gauge-invariant perturbations of a Schwarz- 
schild black hole Q^„^ ■ Overall, we find that both methods 
of calculating the linear momentum loss agree excellently and 
we are thus able to obtain accurate recoil measurements with 
error bars of 5 km/s for the anti-aligned spin binaries and of 
8 km/s in the aligned cases. 
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TABLE IV: Final and radiated angular momenta and masses, computed from the gauge-invariant waveforms. Shown is also the radiated spin 
and mass relative to their initial values, which are listed in Tab.|l] 





0,1 1 12 






r / r 


A/fin 


A/rad 


A/rad/Ai"A 


rO 


-1.00 


0.6244 


0.2008 


0.2434 


0.9536 


0.0320 


0.0325 


rl 


-0.75 


0.6391 


0.2222 


0.2580 


0.9507 


0.0348 


0.0353 


r2 


-0.50 


0.6530 


0.2449 


0.2727 


0.9482 


0.0374 


0.0380 


r3 


-0.25 


0.6676 


0.2670 


0.2857 


0.9461 


0.0396 


0.0402 


r4 


0.00 


0.6827 


0.2886 


0.2971 


0.9439 


0.0420 


0.0426 


r5 


0.25 


0.6966 


0.3106 


0.3084 


0.9412 


0.0450 


0.0456 


r6 


0.50 


0.7075 


0.3363 


0.3222 


0.9376 


0.0488 


0.0495 


r7 


0.75 


0.7181 


0.3626 


0.3355 


0.9344 


0.0523 


0.0530 


r8 


1.00 


0.7292 


0.3878 


0.3471 


0.9315 


0.0557 


0.0564 



Such a good agreement, however, is attainable only if the 
initial transient in the waveform is properly taken into ac- 
count. The transient is produced by the use of initial data 
not containing the net linear momentum the system has accu- 
mulated since inspiralling from infinite separation. We dis- 
cuss the importance of choosing the coiTect vector integration 
constant when calculating the radiated linear momentum and 
describe an unambiguous method for doing so. 

We remark that a proper choice of this constant is essen- 
tial not only because it influences the final recoil velocity with 
differences of 10% and more, but also because it allows for 
a systematic interpretation of the results. Without it, in fact, 
the correct functional dependence of the final recoil velocity 
on the spin ratio is irremediably lost and a comparison with 
the PN prediction impossible. Last but not least, a proper in- 
tegration constant can result in a significant saving of com- 
putational time, allowing simulations to start at much smaller 
initial separations without sacrificing accuracy. 

Tests. In order to show the accuracy of our results, we 
demonstrate that both the Zerilli-Moncrief gauge invariant 
waveforms and the Einstein tensor converge with an order be- 
tween three and four, which is the expected convergence be- 
havior of our numerical methods. 

Furthermore, because the Newman-Penrose scalar 
serves as a measure for the radiation content of the spacetime 
in appropriately chosen frames and at sufficiently large dis- 
tances from the source, we show that the peeling property is 
indeed well satisfied in our numerical simulations. In partic- 
ular, we demonstrate that both the gravitational wave infor- 
mation and the gauge-wave information '^■^ satisfy the ex- 
pected scaling with radius. Similarly, we also show that, as 
expected, the gauge invariant quantity Q ^2 does not vary with 
radius. 

Finally, we investigate those systematic effects that may 
influence our gravitational-wave measurements. In particu- 
lar, we study the effects that the choice of the extraction ra- 
dius has on the final kick velocity and find little influence for 
^ 30 Af. Based on this, we choose = 50 M as the 
fiducial extraction radius in this paper. Furthermore, to ex- 
clude that the effects of the eccentricity in our initial data are 
significant for this paper, we artificially increase or reduce the 
eccentricity of the initial data by comparatively large amounts. 
Also in this case we find that the differences in the recoil ve- 
locities are below the estimated error-bars. Altogether, the set 



of tests carried out gives us confidence that our waveforms 
and recoil velocities are both correct and accurate. 

Results. Using the mathematical and numerical setup as de- 
scribed and tested above, we have investigate the dependence 
of the recoil velocity on the initial data parameters and most 
notably on the spin ratio ai/a2. As expected, a larger asym- 
metry in the initial conditions causes a larger recoil, with a ve- 
locity of about 262 km /s for a binary of equal and anti-aligned 
spins, and a numerically computed recoil of 10^^ km/s for a 
binary of equal and aligned spins. 

Using such accurate measurements, we have then studied 
the functional dependence of the recoil velocity on the spin 
ratio finding that a quadratic behavior reproduces very well 
the numerical results and corrects the post-Newtonian predic- 
tion of a linear dependence. We summarize this behavior in a 
phenomenological expression that can be readily employed in 
astrophysical studies on the evolution of binary black holes in 
massive galaxies. 

With a straightforward extrapolation of the quadratic de- 
pendence to the maximal spinning case ai = —02 = 1 we 
obtain 448 ± 5 km/s as the maximal possible recoil veloc- 
ity attainable from a binary system of equal-mass black holes 
with spins aligned to the orbital angular momentum. This re- 
coil velocity is in very good agreement with our previous esti- 
mate made in Ref. IIQI with a smaller sequence and in equally 
good agreement with the results reported in Ref. . 

As mentioned above, the inclusion of the integration con- 
stant has been essential to obtain physically consistent results. 
At the same time, its investigation has allowed to highlight 
some important features of the evolution of the recoil veloc- 
ity in vector space. Most importantly, it has shown that even 
when all non-spherical modes up to £ ~ 7 are taken into ac- 
count, the recoil is dominated by lower mode contributions, 
especially i ^ 2, m = —2, 1, 2 and £ ~ 3,m = —3. The in- 
terplay of these contributions in vector space and during ring- 
down is what is responsible for the rich features observed in 
the final evolution of the recoil velocity. 

Finally, we provide accurate measurements of the radiated 
energy and angular momentum. These measurements reveal a 
clear linear dependence on the spin ratio a^jax, and we derive 
phenomenological expressions for the relative losses of angu- 
lar momentum and the relative spin-up of the final black hole. 
These relations can be easily used in A^-body simulations if 
the interaction of binary black holes is to be taken into ac- 
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count, and when an accurate estimate on the final distribution 
of black hole spins is important. 
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TABLE V: Integrated convergence rates of the Zerilli-Moncrief 
gauge-invariant variables providing the dominant contribution in the 
kick-velocity measurements. As the numbers indicate, we achieve at 
least third order convergence both in amplitude and phase. A time- 
shift as given by Eqs. l lAU -l lA3t was made on the raw data to remove 
the near cancellation of the lowest-order error terms. 
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APPENDIX B: DETAILS ON THE EXTRACTION OF 1'4 



APPENDLX A: ON THE CONVERGENCE TESTS 

The effects of the initial transient modes can last for dif- 
ferent amounts of time for the different resolutions. A com- 
parison of the (5^2 waveforms between the three resolutions 
confirms this shift in time - the waveform maxima are seen 
at slightly different times for the different resolutions. We 
attempt to undo this effect by manually shifting the time- 
coordinate of the medium and high resolution runs 



t^t + St. 



(Al) 



The value of St is set for the medium and high resolution runs 
independently, using the minimization condition 



d 



d{St) 



170 



150 



\Q{t ^ t + 6t) - Q^higi,\^dt = 0. 



(A2) 



This effectively means aligning in time the peak amplitude of 
the three runs, at <: « 160 M. Solving Eq. (IA2b numerically 
for the Q22 waveforms gives 



StoA 



0.4756 and (5to.oi8 = 0.1078. 



(A3) 



Applying the time-shifting condition Eq. (lAll l to the coarse 
and medium resolution data, and inserting the result into 
Eqs. (I20li-(l2ni gives convergence rates that are consistent with 
the theoretical expectations. 

In Table |V] we report the convergence rates as calculated 
from Eq. |20] for the time interval < u < 190 {u is the 
retarded time as defined in Sec. Ill Eb which excludes the initial 
burst but contains the rest of the waveform. We see close to 
fourth-order convergence for the £ = 2 modes Q22 and Qji- 
The I — m = 3 mode Q^^, on the other hand, shows second 
order convergence in phase, which is most hkely related to 
the fact that the magnitude of this mode is the same size as the 
finite difference error in Q22 and is a factor of 40 smaller than 
the magnitude of Q22 itself. 

The final kick- velocity magnitude in units of km/s is 



|w|kick = 263.49, 259.75, and 261.00 



(A4) 



for the medium, high and very-high resolutions. This gives 
pd'J^lkick) — 2.98 which can be inserted into Eq. d^TTi to obtain 
a calculated convergence rate of 4.32. 



The numerical solution of Eqs. ( |28] l involves first an in- 
terpolation of ^'4 as calculated according to Eqs. (l24l i from 
its values on the Cartesian grid to those onto the extraction 
sphere by using fourth-order Lagrange interpolants. Because 
of the symmetry across the z = plane the interpolation is 
effectively done on the upper hemisphere only, thus using a 
spherical coordinate system with 6*, G [0, 7r/2] x [0, 2tt] and 
applying cell-centered discretization along the ^-direction to 
avoid the coordinate singularities at the poles on the sphere. 

The angular resolution is chosen so that the spacings A9 
and Ac/) are equal and of the same order as the corresponding 
Cartesian spacings of the refinement level in which the largest 
extraction 2-sphere is located. As an example, for the fiducial 
finest resolution of h = 0.024 M, the largest extraction radius 
is at Tj, = 60 AI and in a region covered by the second re- 
finement level with spacing A'^1'^2 = 1-536 M. To obtain an 
equivalent spacing on the 2-sphere, we solve for A6 and A(f> 
such that 



,A0 



0.024 
2 



1.536 M. 



(Bl) 



The resulting number of grid points is Ne = 56 along the 
6'-direction and = 224 along the (/)-direction. 

After interpolation onto the extraction sphere, we first cal- 
culate the time integral of 5*4 1 52 and afterwards, the surface 
integral of the absolute square of the former according to 
Eqs. ( |28] ). These integrals are both computed using fourth- 
order schemes. In particular, for the surface integral, we use 
Simpson's rule in the form 



dxf{x) 



Aa 



17 59 43 49 



49 43 59 17 ■ 



, (B2) 



where (fk) is the sum over all fk with 3 < fc < — 3. 
The integral over ddd(j> is obtained by computing the tensor 
product of the RHS of Eqs. dMI i. i.e.. 
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FIG. 13: Left panel: Evidence that the conditions for the Peeling theorem are met also for ^'3, which scales as r^^ when extracted at isotropic 
radii r j, = 30 M, 40 M, 50 A/, and 60 M. This figure should be compared to the corresponding Fig. |4] Right panel: The same as the left 
panel but for the gauge-invariant quantity Q^2' which is shown to be constant when extracted at isotropic radii r^, = 30 M, 40 M, 50 M, and 
60 M. 



where the c;, Cj are the coefficients in the RHS of Eqs. 

The time integral of Eqs. ( |28] | is generically calculated by 
using the fourth-order Simpson's rule in such a way that the 
integral for the time step k uses only past time steps i with 
< i < fc. Care is required for the very first time steps, for 
which we have less than 7 evaluations of the integrand. In this 
case, we use the 2nd-order accurate trapezoid rule if 1,3, 
or 5 



dx / (a 



Aa 



+ (./fe) + -j^fN 



(B4) 



or the fourth-order accurate Simpson's rule 



dxf{x) 



Ax 



3/0 + 3/1 



2 4 1 

+ (3/2*; + 3/2^+1) + 3/jv 



, (B5) 



if iV — 2, 4 or 6. For > 7 we simply use Simpson's rule in 
the form ( IB2I ). It should be noted that the use of a higher-order 
time integration scheme improves the overall accuracy in the 
calculation of the final recoil velocity by more than a factor of 
10. 



APPENDIX C: A COMPARISON OF WAVE-EXTRACTION 
METHODS 

In Fig. m we have shown that vI/4 as extracted at different 
radii correctly scales with the 1/r falloff as predicted by the 



peeling theorem. Here, we also check if all other components 
of the Weyl tensor exhibit the correct r^^^'i'n = const, scal- 
ing. 

The left panel of Fig.[T3]indeed shows that the scaling prop- 
erty of all behave as expected. In the course of the same 
analysis, it is also worth looking at the waveforms as calcu- 
lated by using the gauge-invariant formalism. In particular, 
we focus on the real part of the £ = 2,m = 2 even parity 
wave mode Q22 and check for the correct scaling for the dif- 
ferent extraction radii. The right panel of Fig. [13] shows that 
Q22 is constant for all extraction radii as expected. 

As a final remark, we will also compare the and /ix as 
calculated by using the odd and even master functions in the 
gauge-invariant formalism according to Eq. ( [3Q] l and the spin- 
weighted spherical harmonic amplitudes of the Weyl compo- 
nent decomposed on the extraction spheres. Using these 
amplitudes, the metric perturbations /i+,/ix recovered by a 
double time integral of Eq. 



— i/ix = lim > 

r — ^00 ^ — ^ 



dt' 



dt"^' 



(CI) 



The numerical integration of Eq. dClt requires knowledge of 
an integration constant for the calculation of the second in- 
tegral to eliminate the linear offset. This constant is deter- 
mined by searching for minima in the ^'4™ mode and averag- 
ing over them. The resulting value is used as the integration 
constant. In both cases, we only consider the dominant con- 
tribution from mode £ = 2 , m = 2. 
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FIG. 14: Comparison of the two polarization amplitudes h+ (up- 
per graph) and h-^ (lower graph) as computed with 3*4 (continuous 
black line) or with the gauge invariant quantities (dashed red 
line). Note the two polarizations are computed using the lowest (and 
dominant) multipole 1 = 2, m = 2 and are extracted at = 50 M. 



evolutions begin from fairly close separations, comprising at 
most the last 2-3 orbits. As such, parameters for quasi-circular 
orbits determined by the effective potential method, give only 
approximations to the true orbital parameters for black holes 
that have spiraled in from infinity, and it is known that the 
method produces a non-trivial residual eccentricity for initial 
data at close separation. This eccentricity can have signif- 
icant effects on the orbital trajectories before merger, and a 
potential influence on the calculated recoil. To test this we 
have evolved two modified rO models, one in which the ini- 
tial linear momenta of the black holes is 3% larger than that 
specified in Table U and another in which the linear momenta 
are 3% smaller. The modified momenta have the effect of 
changing the orbital energy of the bodies from the minima 
determined by the effective potential method, introducing an 
additional eccentricity to the evolution. The resulting black 
hole trajectories and kick determinations are shown respec- 
tively in Fig. [15] We see that although the level of applied 
eccentricity is large, and in fact much larger than the expected 
eccentricity due to the intrinsic inaccuracy of the effective po- 
tential method, it modifies the recoil by only about lOkm/s, 
that is, 4%. Further, in both the high and low energy cases, 
the recoil is increased over the fiducial rO case, suggesting that 
increased eccentricity generically leads to a slightly larger re- 
coil. 



APPENDIX D: ON THE INFLUENCE OF ORBITAL 
ECCENTRICITY 

Another source of potential error in calculating a "physical" 
kick comes from the choice of initial data parameters. Our 
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